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We present theoretical approaches to high energy nuclear collisions in detail putting a 
. . . special emphasis on technical aspects of numerical simulations. Models include relativistic 

04 \ hydrodynamics, Monte-Carlo implementation of fcT-factorization formula, jet quenching in 

expanding fluids, a hadronic transport model and the Vlasov equation for colored particles. 

o 

^^ ! §1. Introduction 

One of the primary purposes of heavy ion physics at ultrarelativistic energies is to 

explore properties of strongly interacting matter at high temperature and densities, 

^^ ' namely, the quark gluon plasma (QGP).'^ ^'tU To this purpose, collider experiments 

of high energy pp, dA and AA collisions have been performed at both Relativistic 

Heavy Ion Collider (RHIC) in Brookhaven National Laboratory (BNL) and Large 

Hadron Collider (LHC) in European Organization for Nuclear Research (CERN). It 

is of particular importance to comprehensively understand space-time evolution of 

^ \ matter created in high energy heavy ion collisions so that one can extract information 

Ch ' on the detailed properties of matter from experimental data. 

High energy nuclear collisions contain rich physics and exhibit many aspects of 

dynamics according to relevant energy and time scales. Collisions of two energetic 

'^ ' nuclei can be viewed as those of highly coherent dense gluons. Universal behaviors 

of hadrons and nuclei at the high energy limit are called the color glass condensate 

'sj" \ (CGC).^' Just after the collision, longitudinal color electric and magnetic fields be- 

''^ ' tween the two passing nuclei, which is called color flux tubes,'^ are produced from 

rn . CGC initial conditions. Subsequent evolution of these flux tubes is called "glasma".!^ 

Long range rapidity correlations in the glasnicP'''^ ™ provide a natural explanation 

for the ridge structure observed at RHICP^''^ and LHC.^ It was pointed out that 

instabilities of the gluon fields could play a significant role in the process of thermal- 

, ization-ISJ-nS-HB 

K< ■ A model based on fcr-factorization iormal^^^^^^^^^ or clas- 

Cd . sical Yang- Mills approachP3'l2Sli reproduces the multiplicity distribution for charged 

hadrons at RHIC and LHC. However, initial transverse energy per particle is large 

compared with the experimental data.'^''^''22i Qlasma evolution based on a 2-1-1 

dimensional classical Yang-Mills simulation does not account for elliptic flow data 

at RHIC.^^i* These facts suggest the necessity of inclusion of further evolution with 

much stronger interaction, e.g., hydrodynamical evolution. Indeed, a nearly perfect 

fluid pictur^''33''E3''321''l33 turns out to explain large elliptic flo-KiSl observed at 

RHIC'^^'"^^''''^ and LHC,'3SJS'linii which leads to establishment of a new paradigm 

typeset using PTPTeX.cIs (Ver.0.9) 
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"strongly coupled QGP (sQGP)"P''3'''33'H3 These hydro dynamic simulations re- 
quire very short {^ 1 fm/c) thermalization time. 

Due to expansion and cooling down, a QGP fluid becomes eventually a hadronic 
gas at a late stage whose evolution can be described by hadronic transport mod- 
g^gp3|i,|3l]i„44),45),46) j^galistic gradual freeze-out both chemically and kinetically can 
be naturally treated in these models. It is claimed that a QGP fluid picture and a 
hadron gas picture are demanded to understand pT speactra and differential elliptic 
flow parameter for identified hadrons simultaneously;^^' For the hadronization pro- 
cesses, it has been claimed that quark recombination/coalescence is important for 
intermediate transverse momentum region.SZJ 

On the other hand, high momentum jets are also created in the collision at 
collider energies. These jets have to traverse the matter in heavy ion collisions 
and, therefore, interaction with the expanding matter should be treated accordingly. 
During traveling through the medium, jets interact with soft matter and lose their 
energies (jet quenching). Therefore, high transverse momentum hadrons are good 
probes of the bulk matter .SSJ 'SSI .l5QJi,|5]j Suppression of high transverse momentum 
hadrons was observed at RHIO^''^ and LHC.'^ Disappearance of the away-side 
peak in azimuthal correlation functions for high transverse momentum hadrons is 
also observed in central Au+Au collisions at RHIC. ■ High pT suppression in heavy 
ion collisions turns out to be attributed to the final state interaction because of the 
Cronin enhancement and the existence of back-to-back correlation in d+Au collisions 
at midrapidity at RHIC.'^ Theoretical approaches to describe such jet quenching 
need space-time evolution of parton density through a trajectory of a jet. Hydro- 
dynamical simulation provides such a parton density. A hydro + jet model was 
proposed in Refs. [57j) . l58]) later followed by Refs. [59]) . [60|1 . [6T]) incorporating more 
realistic models for parton energy loss in a dense medium. 

Recent experimental data on higher order anisotropic flow^'ISS'lSIli c^ll atten- 
tion to the importance of initial state fluctuations in heavy ion collisions. It is 
known that event-by-event fluctuationgSS lead to the higher order anisotropic flow 
such as triangular flow quantified by third harmonic component of azimuthal angle 
distributions, and importance of event-by-event simulations with fluctuating initial 
conditions has been realized. '^^''^^''^^''^^''^^''^''^ The triangular flow contributes 
also most of the ridge observed in the two-particle azimuthal correlation. So far only 
fluctuations from configuration of nucleons inside a colliding nucleus are included in 
the hydrodynamical simulations. Recently fluctuations from particle production it- 
self are studied^ by including negative binomial distribution which was obtained by 
the glasma.l^ It was found that multiplicity fluctuations in pp and pA collisions at 
midrapidity exhibit Koba-Nielsen-Olsen (KNO) scaling.^^ Fluctuations of particle 
production increase higher order anisotropy such as triangularity by about 50%. '^ 
DIPSY event generator which includes the fluctuations arising from dipole evolution 
also predicts larger higher order anisotropy.l^ Viscous hydrodynamic simulations 
with fluctuating Glasma initial conditions have been performed recently.'^ Monte- 
Carlo Glauber type model in which effects of nucleon-nucleon correlationaZSJ-lIO ^j-g 
included was proposed to construct more realistic nuclear configurations for initial 
conditions of subsequent evolution. 
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Fig. 1. Dynamical modeling of relativistic heavy ion collisions in view from proper time and energy 
scale. 



As described above, it is needed to incorporate all such different physics consis- 
tently to have unified and better understanding of the space-time evolution of the 
system created in high energy nuclear collisions. Figure [T] shows several important 
aspects of dynamics of relativistic heavy ion collisions according to time and energy 
scales. Experimental observables reflect all the history of evolution of matter staring 
from initial colliding nuclei to final free-streaming hadrons. A first attempt to an in- 
tegrated approach to the heavy ion collision as a whole was done in Ref. [78jl in which 
full three dimensional ideal hydrodynamic simulations with initial conditions taken 
from a CGC based model were performed and parton energy loss was simulated in 
these expanding fluids. In this review, we present technical and numerical aspects 
of these important modules for high energy heavy ion collisions; CGC, relativistic 
hydrodynamics, parton energy loss, hadronic transport models and Vlasov model for 
colored particles. 

This paper is organized as follows. In Sec. [21 we explain in detail numerical 
aspects of relativistic ideal hydrodynamics in heavy ion collisions. In Sec. [3l Monte- 
Carlo implementation of the CGC initial conditions based on the /c^-factorization 
formula is discussed. Energy loss of energetic partons in an expanding QGP fluid 
is briefly discussed in Sec. HI In Sec. [5l a cascade method and cross sections in the 
hadronic transport model JAM are briefly summarized. In Sec. El we discuss how 
to solve the Vlasov equation for colored particles by employing the particle-in-cell 
method. The final section is devoted to summary and conclusion. 



§2. Relativistic Hydrodynamics 



Relativistic hydrodynamics is one of the key dynamical framework to describe 
the space-time evolution of matter created in relativistic heavy ion collisions. Since 
the main goal in the physics of relativistic heavy ion collisions is to understand 
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the properties of matter under (local) equilibrium, one can apply hydrodynamics, 
in which local thermal equilibrium is assumed, to dynamical description of created 
matter in any cases as a bottom-up approach to see whether the hydrodynamic 
description works well. In this section, we briefly overview framework of relativistic 
ideal hydrodynamics and its numerical aspects. 

2.1. Relativistic ideal hydrodynamics 

Relativistic hydrodynamic equations describe conservation laws of energy and 
momentum 

a^r'^'^(x) = 0, (2-1) 

together with the conservation of charges 

d^Ntix) = 0. (2-2) 

Here T^'^ is the energy momentum tensor and N-^ is the i-th conserved current. With 
an assumption of ideal hydrodynamics where all dissipative effects are neglected, one 
can decompose the energy momentum tensor and the conserved currents as 

Tf"" = euf'u'' -Pigf"" -u^'u"), (2-3) 

K = n^u^^ (2-4) 

where e, P, rii and u^ = 7(1, i^) = -^=^{l,v) are energy density, pressure, i-th 
conserved charge density and four flow velocity, respectively. Minkowski metric in 
this paper is defined as (7^^ = diag(l, —1, —1, —1). Instead of Eqs. (|2-ip and (|2-2p . 
the following expression of the balance equations might be also convenient when 
compared with non-relativistic equations: 

^^E + V-iE + P)v = 0, (2-5) 

^M'' + V • M'^v = -V'^P, (2-6) 

^N, + V- N,v = 0, (2-7) 

at 

where 

E = {e + P)j^ - P, (2-8) 

M = {e + P)j'^v, (2-9) 

Ni = na- (2-10) 

In ideal hydrodynamic framework, the equation of state plays an important 
role. First, the hydrodynamic equations ()2-ip and ()2-2p are not closed as a system of 
partial differential equations: the number of unknowns is 6 (energy density, pressure, 
charge density and three components of flow velocity) in the case of one conserved 
charge while the number of equations is 5. So the system is closed when a relation 
among unknowns, e.g., P = P{e,n) is specified. Second, equations (|2-ip and (|2-2p 



just describe conservation laws and are the so-called balance equations. Under an 
assumption of local thermal equilibrium, one can utilize the equation of state P = 
P{e, n) which only reflects the microscopic dynamics. Since the collective flow is 
generated by pressure gradient, it is sensitive to the equation of state and the degree 
of kinetic equilibrium. This is one of the reasons why the collective flow has been 
focused in relativistic heavy ion collisions. 

The equation of state can be taken from results from the first principle calcula- 
tions of QCD thermodynamics, namely, lattice QCD.'^'I^''^''^ A few comments 
regarding this are in order here. 

1. In the current status of the lattice QCD results, the equation of state is not so 
reliable in the low temperature region (T^lOO MeV). Then one can connect lat- 
tice QCD results at high temperature with results from the hadronic resonance 
gas model at low temperature. ^^''^^* Practically, this is also demanded from a 
point of view of freezeout at which all hydrodynamic variables are switched to 
a particle picture employing the Cooper-Frye formula.!^ Energy, momentum 
and charges are conserved in this formula only when the hadronic resonance 
gas picture is valid. 

2. Monte-Carlo calculations of lattice QCD at finite baryonic chemical potential^ 
suffer from a severe issue due to the so-called sign problem. The equation of 
state at finite baryon density, however, would be demanded in lower collision 
energies or in forward/backward rapidity regions. 

3. Even if results will become reliable in the low temperature region and/or finite 
baryon density region in lattice QCD, there is an issue on chemical freezeout 
since all thermodynamic variables are obtained for thermally as well as chem- 
ically equilibrated states. In the actual hadronic matter created in relativistic 
heavy ion collisions, chemical composition of hadrons is almost frozen during 
expansion according to statistical model analyses. Thus, each hadron acquires 
its chemical potential associated with the approximated conserved number of 
the hadron below chemical freezeout temperature. Again, the hadronic res- 
onance gas model with finite chemical potential^''^''^''^''^ is needed to 
describe the space-time evolution of hadronic matter in hydrodynamics. Since 
the matter is already diluted due to expansion, one can instead use the hadronic 
cascade model for a better dynamical description of hadronic matter. This will 
be discussed in Sec. [H 

2.2. Numerical aspects of relativistic ideal hydrodynamics 

In this subsection, we review a numerical scheme to solve relativistic hydro- 
dynamic equations for perfect fiuids. We first discretize relativistic hydrodynamic 
equations in Cartesian coordinate to solve numerically in Sec. 12.2.11 We also discuss 
the conventional way to treat multi-dimensional problem, namely, the operator split- 
ting method. After that, we introduce the piecewise parabolic method (PPM),'^ 
which is known as a robust algorithm against strong shock waves, to relativistic 
hydrodynamic equations in Sec. 12.2.21 In practice, one needs to convert a set of 
numerical solutions to physical quantities such as thermodynamic variables and flow 
velocity. The procedure to obtain them is explained in Sec. 12.2.51 Finally, we discuss 
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hydrodynamic equations in relativistic coordinate (proper time r and space-time 
rapidity rjs) in Sec. I2.2.4I 

2.2.1. Discretization and the operator splitting method 

In the Cartesian coordinate, hydrodynamic equations ()2-ip and (|2-2p for a rela- 
tivistic perfect fluid with one conserved charge [e.g., baryon charge) can be written 
as 



dt 



where 
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(2-11) 



(2-12) 



These equations can be summarized as a continuity equation in the following form 



dtUjit,x)+ Yl d^Fu[Ujit,x)]=0, (J =!,•••, 5). 



(2-13) 



i=x,y,z 



For example, F^i = UiVx + P, Fyi = UiVy, Fx2 = U2VX and so on. 



One discretizes Eq. (|2-13p and write in a general form 

~Ax 



[U. 



in+l 
Hjk 



[U. 



Jiijk 



l=x,y,z 



Gij [[UjY^jk] 



or more simply. 



^r 



Tjn 



Ax 



Gi [U^jk 



(2-14) 



(2-15) 



Here n is a time step and z, j and k are fluid cell indices in x, y and z direction, re- 
spectively. At and Ax are mesh sizes in temporal and spatial direction, respectively, 
which should obey C < 1 in relativistic cases where C = At/ Ax is the Courant 
number (Note that we have employed natural unit c = 1). We here assume isotropic 
lattice in three dimensions. One can cope with these kinds of multi-dimensional equa- 
tions by employing the operator splitting method. In this method, the operators are 
split into three sequential one-dimensional spatial steps: 



'Jijk 
^ijk 



US' 



'-^ijk 



'-'ijk 
^ijk 



■^GyPijk], 



At_ 
~Ax 



Gz[Uijk\- 



(2-16) 
(2-17) 
(2-18) 



To avoid numerical errors on spatial anisotropy, the above process is cyclically 
changed every time step. Now the problem in three dimensional space reduces to 
the one in one dimension. Hereafter in this subsection we concentrate our discussion 
on solving hydrodynamic equation in one-dimensional space. 

2.2.2. Piecewise parabolic method 

As a robust numerical algorithm to solve hydrodynamic equations, we review the 
piecewise parabolic method (PPM)^^' in this paper. PPM was employed for the first 
time in the physics of relativistic heavy ion collisions in Ref. [92]) to solve relativistic 
hydrodynamic equations in Cartesian coodinate and later applied to problems in rel- 
ativistic T-r]s coordinate in Ref. [32j). PPM is categorized in the so-called Godunov's 
method. Godunov made use of a Riemann's shock tube problem, approximated a 
solution of interpolating two discontinuous states to a constant one and obtained it 
using conservation equations. PPM is, so to speak, a higher order extension of the 
Godunov's approach. As a consequence, PPM allows one to describe a hydrody- 
namic response to steep profile very efficiently. For details, see Ref. I9ip . For other 
algorithms such as SHASTA and rHLLE, see also Ref. |93]) 

For a given set of discrete values of fields {Uj"} at time step n, an interpolation 
function, Uj{x), can be defined as 

I f^i + l/2 



[/; = —/ U,ix)dx, (2-19) 

^^ •''a;j-i/2 

where Uj{x) is assumed to be continuous and a parabolic function in Xj — Ax/2 = 
Xj_i/2 < X < Xj+i/2 = Xj + Ax/2. Then one can parametrize Uj{x) as 



AUj = Ur,j-Ulj, (2-21 



X — X,_l /9 

AU, + Ue,{l ^^ 



U6,i = 6 



U^-l^iURJ + ULj] 



2-20) 

2-21) 

(2-22) 



where Uj{xj_i/2) = Ulj and Uj{xj^i/2) = Ur,j to be determined as follows. 
As default values at cell boundary, one puts 



ur, = i^(f^r + f^iVi) - Y^(f^;v2 + uf-i), (2-23) 

UL,i = ^(C/jLi + UJ) - ^{U^+1 + U^), (2-24) 



using discretized solutions {U!/} at time step n. These values are obtained as follows. 
Assuming integral of U{x) can be parameterized using a quartic function (here U{x) 
is rather defined globally at least in Xj_^/2 < x < Xj^^/2), 



U{x) = j U{x')dx' (2-25) 

= ax^ + bx^ + cx^ -h dx -F e, (2-26) 
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we calculate the value at the cell boundary x = Xjj^in as 

(2-27) 



TT( ^ ^ 



dx 



^i + l/2 



Only when we solve this problem, one can suppose a;j+i/2 = and U{x = Xj_3/2) = 
without loss of generality. Then, U{xjj^ii2) = d. We obtain the following coupled 
equations, 

U{2Ax) = / U{x)dx = {Uj^i + Uj + Uj+i + Uj+2) ^x, (2-28) 

•^^j-3/2 

U{Ax) = (C/j_i + Uj + Uj+i)Ax, (2-29) 

U{0) = {Uj-i + Uj)Ax, (2-30) 

Ui-Ax) = Uj-iAx, (2-31) 

U{-2Ax) = 0. (2-32) 

Solving these equations with respect to a, • • • and e, we finally obtain Eq. ()2-23p . 

However, the value has to be reset in some cases. The interpolation function 
Uj{x) is imposed to be monotonic in Xj_i/2 < x < Xj_^_i/2- Uj{x) should take its value 
between [//jj and Ul,j- In fact, Uj{x) does not obey this condition either when C/J' 
is a local minimum or maximum or when Uj{x) has an extreme even though [/" is 
between Urj and Ulj- From Eq. (|2-2U|) . this is the case when | AUj \>\ Uqj \. The 
following replacement is made in these cases: 

Ur,j,Ulj^U^ if {Ub,j-U;^){U^-Ul,j)<0, (2-33) 

Ul,j -^ 3C/7 - 2Urj if AUjUe,j > {AUjf, (2-34) 

Urj -^ 3U^ - 2Ul,j if ^UjUej < -{AUjf. (2-35) 

In the case of heavy ion collisions, fluids sometime can be surrounded by vacuum 
in which U^ vanishes. In this case, both Ulj and Urj are set to be zero. When the 
default value calculated using Eq. (|2-20|) becomes negative, it is also set to be zeror^l 

Using Ulj and Uptj determined above, one calculates thermodynamic variables 
{&L,R and Pl,r), sound velocity {cl,r) and flow velocity {vl,r)- We will discuss how 
to obtain these variables from numerical solutions U in the next subsection. Average 
values of interpolation function around the cell interfaces are 

^L,j+i = TT nr. / C/,+i(x)dx, (2-36) 



K,j+l/2 I ^t 7^^.^j/2 

/rXj+1/2 
.. . , ■ ^j + l/2~Kj + l/2\^t 



0; o^i/o Zit ' 



*' This might have been too strict for t/5 to be set since it could be negative. However, if initial 
f/5 is positive, it keeps to be positive in this fluid element. Therefore there is no problem in ordinary 
cases. 



br and bi are signal velocities 



&rj + l/2 



max I 0, 



VL,j+l + CL,j+l ^i+1/2 + Cj+1/2 



1 + VL,j+lCL,j+l 1 + ^^+1/20^+1/2 
;, - ^^^ I n ^RJ - ^RJ ^i+1/2 + 9+1/2 ^^ 

bij+1/2 - mm ( 0, ,_^,__, , T^^r—Z^—^j > 



1 - VRJCRJ 1 + t;j + i/2Cj + i/2 , 



V-i 



i+1/2 



Cj+1/2 



(Uijj + VL. 



J+lJ' 



(Ci?,j+CLj+i). 



(2-38) 
(2-39) 

(2-40) 
(2-41) 



Notice that bijj^i/2 < 0. When a fluid cell is located next to vacuum, e.g., {Ui)j ^ 
and (C/4)j+i = 0, we set br = 1. Please notice that the subscripts R and L represent, 
respectively, the values of right side and left side at each cell, but that the subscripts r 
and I denote, respectively, the values of right side and left side at each cell boundary. 
Inserting Eq. ([^^ into Eqs. (^^ and ([2^ . we obtain 



Ul,j+i = Ul,j+i + 
Urj = Urj 



br,j+l/2'^t 



2Ax 
2Ax 



AU. 



i+i 



AU, 



+ 1 



1 



26. 



r,j+l/2 



At 



26,, 



';,i+i/2 



3Ax 
At 



U> 



6,i+l 



3Ax 



Ua 



(2-42) 
(2-43) 



By using these average values above, we solve the Riemann problem at each cell 
boundary. The solution becomes complicated in general. Thus, in the Godunov- 
type algorithm, one approximates the solution to a constant value Uir which fulfills 
the conservation law. One rewrites the hydrodynamic equations in their integral 
form: 

[U{x, tn + At/2) - U{x, tn)]dx 

"-3 



tr,+At/2 



The integration can be easily done by assuming U{xj] 
tn + At/2. The result becomes 

^* -[F{U,,)-F{Ur,)\ 



[F{U{xj+i/2,t))-F{U{xj,t))] dt (2-44) 

Urj = const, in i„ < t < 

At 



One also integrates the hydrodynamic equations over Xj+1/2 < x < Xj+i 

Xj + l 

[U{x, tn + At/2) - U{x, tn)]dx 

^J + 1/2 

t„+At/2 

[F(?7(x,+i,t)) -F([/(x,+i/2,t))] dt 



(2-45) 



and obtains 



{Ulr - UL,Hl)br,j+l/2^ = " [F{Ul,j+i) " F{Ulr)] ^. 



(2-46) 



(2-47) 
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From Eqs. (12^151) and (l2^T7l) . one obtains Uir and F{Uir) 

F{Urj) - F{Ul,j+i) - bij+i/2UR,j + brj+i/2UL,j+l 



F{Ulr) 



K 



i+i/2 



F{Ur,) 



Vi+1/2 ~ ''/,j+l/2 



h,j+l/2P{UL,j+l) + &rJ+l/2^/,j+l/2 



(2-48) 



(Ulj+i - Urj] 



\,j+l/2 — "l,j+l/2 



Finally, the solution at the next time step leads to 



u: 



n+l 









-1/2 



F, 



J-l/2) 



(2-49) 



(2-50) 



The key difference between the conventional Godunov scheme and the present method 
is using the average values Urj and Ul,j+i instead of Uj and C/j-+i to gain accuracy of 
numerical solutions since numerical fluxes using Uj or f7j+i are often overestimated, 
in particular, in the case of steep proflles. 

2.2.3. Thermodynamic variables and flow velocity 

We transform from thermodynamic variables and flow velocity to variables to 
solve the relativistic hydrodynamic equations numerically in Eq. ()2-12p . Hence we 
need to transform back to thermodynamic variables from numerical solutions Uj. 
From Eq. (|242]) . 

\U\ 



V 



riB 
U 



U4 + P{e,nB)' 
Ua-U -v, 



iUi,U2,Us) 



(2-51) 

(2-52) 

(2-53) 
(2-54) 



For a given equation of state, P = P{e, ns), and numerical solutions, Uj, Eq. (|2-5ip 
with Eqs. ()2-52p and ()2-53p becomes a non-linear equation with respect to v. There- 
fore, one has to solve it numerically using, e.g., the bi-section method. Once the 
solution V is obtained from Eq. (|2-5ip . it is easy to obtain energy density and baryon 
density from Eqs. (|2-52p and (|2-53p and, consequently, pressure from the equation of 
state. 

2.2.4. Hydrodynamic equations in relativistic coordinate 

It is more appropriate to write down relativistic hydrodynamic equations in an 
expanding coordinate in the case of relativistic heavy ion collisions: 
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0,(2-55) 
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where 



/Ui\ 
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T^^{e + P)vy 

T^\e + P)v, 

T^^{e + P)-TP 



(2-56) 



where proper time r = \/t'^ — z"^, space-time rapidity rjs = (l/2)ln[(t + z)/{t — z)] 
and V = {dx^dy^dri^/r). Flow velocities and fluid rapidity in this coordinate are, 
respectively, 



cosh Yi 



f 



-"ris 



cosh(17 - r]s) 

cosh Yf 
cosh(yj — "qs) ^' 
tanh(y/- - r?s), 
1 



1 



V — V — V 



(2-57) 

(2-58) 
(2-59) 
(2-60) 



where ^7 = | In [(1 + f^)/(l - v^)]. 

Equation (|2-55p is quite similar to Eq. (|2-11|) except for existence of the last term 
in the left hand side which is a source term due to expanding coordinate. Therefore, 
one can utilize the same PPM algorithm by adding source terms in the cycle of 
operator splitting. 

One can obtain thermodynamic variables and flow velocities from the following 
relations: 



U\ It 



U^/t +p{e,nBy 
Ua \U\ , ^ , 



UB 



T 
U5 
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(2-61) 

(2-62) 
(2-63) 



Inserting Eqs. (j2-62p and (j2-63p into Eq. (j2-6ip . we first solve an implicit equation for 
I V I numerically and then obtain e, hb and P = P{e, ub)- From numerical solutions 
Uj, we also obtain each component of v: 



1 



V 



-,iUl,U2,U3). 



Ui + rP' 

Velocities in the Cartesian coordinate are obtained from v 

cosh.{Yf — rjs 



cosh Yf 
cosh (17 -r]s) . 



-Vrr-, 



cosh Y) 



f 



yy, 



(2-64) 

(2-65) 
(2-66) 
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Vrj^ cosh rjs + sinh rjs 



tanhy/, (2-67) 

cosn r]s + Vj^^ smn ry^ 

cosh Yr 

cosh(y/- - r/s) 



where Y/ = tanh~ li^^ + r/^ . 

§3. Initial Conditions 

We need to specify initial conditions for hydrodynamic simulations. In principle, 
we must understand the initial particle production and subsequent non-equilibrium 
evolution of the system toward thermal state to obtain the initial condition. However, 
at the moment, we do not have complete understanding at early stages of high energy 
nuclear collisions. Here, we make a simple assumption that produced particles after 
the collision of two nuclei can be used to obtain initial entropy distribution. We 
will come back to this issue on isotropization at an early stage of collisions from a 
viewpoint of non-Abelian plasma instability in Sec. [6l 

The /cy-factorization formulation is widely used to compute the inclusive cross 
section for produced gluons^''^'^''^'^ in hadronic collisions. In order to study 
gluon production in nucleus-nucleus collision, we shall use the Monte-Carlo imple- 
mentation of /cT-factorization formulation or Glauber model (MC-Glauber)24l)''2"''i2Sl' 
in which fluctuations of the position of nucleons inside a nucleus are taken into 
account. Thus we can study nucleus-nucleus collisions on an event-by-event basis. 

We first sample the positions of nucleons inside a nucleus according to a nuclear 
density distribution {e.g., Woods-Saxon function) for two colliding nuclei, and shift 
them by a randomly-chosen impact parameter h with probability bdh for an event. 
A nucleon-nucleon collision takes place if their distance d in the transverse plane 
orthogonal to the beam axis fulfills the condition 

(3-1) 

where Uin denotes the inelastic nucleon-nucleon cross section. Incident energy de- 
pendent total pp cross section is parameterized by Particle Data Group .1^ Elastic 
cross section is computed using PYTHIA parametrization.'^''^ The following val- 
ues are obtained; CTin = 39.53,41.94 and 61.36 mb at ^ = 130,200 and 2760 GeV, 
respectively. In this way, we obtain the number of binary collisions and that of 
participants for each event. It should be noticed that the standard Woods-Saxon 
parameters shown in, e.g., Ref. 1100]) cannot be directly used to distribute nucleons 
inside a nucleus because of the finite interaction range in our approach. We need to 
modify nuclear density parameters so t hat a convolution of nucleon profiles leads to 
the measured Woods-Saxon profile.'^^ll 

Next, we compute particle production at each grid in the transverse plane. In the 
MC-Glauber approach, we assume that the initial entropy profile in the transverse 
plane is proportional to a linear combination of the number density of participants 




and that of binary collisions: 

dS 



so{x±) 



Todxdydrjs 
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= — ( ^^— /3part(r±) + apco\\{ri_) ) , (3-2) 



where tq = 0.6 fm/c is a typical initial time for the hydrodynamical simulation. Pa- 
rameters C = 19.8 and a = 0.14 have been fixed through comparison with the cen- 
trality dependence of multiplicity data in Au+Au collisions at RHICpSS by pure hy- 
drodynamic calculations with temperature Tjec = 100 MeV.'^^ At the LHC energy, 
C = AlA and a = 0.08 are chosen^"^' so that we reproduce the ALICE data on cen- 
trality dependence of multiplicity in Pb+Pb collisions at t/saF/v = 2.76 TeV.123J''ll25J 
The participant density Ppaxtif^) at each grid point is the sum of participants 
density pA^fi.) from nucleus A and psi'^A.) from nucleus S, which are computed 
by counting the number of wounded nucleons Na,w and Nb,w for nucleus A and B 
within a tube extending in the beam direction with the radius r = \J oi^j-k: 

/Opartl^-x) = /9a(^±) + /Ob(^±) = —^ (3-3) 

Similarly, the number of binary collision density at each grid is obtained by counting 
the number of binary collision A'coii with the area din, where the transverse position 
of binary collision is assumed to be the average transverse coordinate between two 
colliding nucleons 

Pcoll(r-x) = ^^ . (3-4) 

fin 

which may be also obtained by the expression PA(^±)p_B('"±)o"in- 

In the Monte-Carlo KLN (MC-KLN) model,!^ the number distribution of gluon 
production at each transverse grid is given by the /c^-factorization formula^*^' 

dNg 4iVe /d^ .^2^^ 



d^r^dy N^ -1 J p\ j 4 

X ^a{xu {p± + fc±)V4) <pBix2, (p±-fc±)V4) , (3-5) 

with A^c = 3 the number of colors. Here, p± and y denote the transverse momentum 
and the rapidity of the produced gluons, respectively. The light-cone momentum 
fractions of the colliding gluon ladders are then given by xi_2 = p± exp{±y) / ^/Fmn , 
where ^/snn denotes the center of mass energy. Running coupling as{Q'^) is evalu- 
ated at the scale Q^ = max((p_L — k±)^/A, {p± + k±)'^ /4). An overall normalization 
factor K is so chosen that the multiplicity data in Au-|-Au collisions at RHIC are fitted 
in most central collisions. In the MC-KLN model, saturation momentum is param- 
eterized by assuming that the saturation momentum square is 2 GeV^ at x = 0.01 
in Au+Au cohisions at 6 = fm at RHIC where ppart = 3.06 fm~^!23 

QlA.:r,) = 2GeV'J^J'^)\ (3.6) 

1.53 tm \ X J 
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where A is a free parameter which is expected to have the range of 0.2 < A < 0.3 
from Hadron Electron Ring Accelerator (HERA) global analysis for x < o.Ol.liSa in 
MC-KLN, we assume the gluon distribution function as 



asiQlj^)max{Ql^,kj 



, 2. ^ -■- ^s,A 



4>Aix,kj_;rj_) r^ — o s 77^ TTT ' i^''^) 



We assume that initial conditions of hydrodynamical simulations are obtained 
by identifying the gluons' momentum rapidity y with space-time rapidity rjg 

so(a;±) oc — -— (3-8) 

Note here that recently gluon distribution function obtained from numer- 
ical results of the running coupling Balitsky-Kovchegov (rcBK) evolution equa- 

^^Qj JT071l . [T08ll. [l09ll 

dj\ [T X] 

"■ ' ' - I d^n K{r,ri,r2)[M{ri,y) + M{r2,y) - M{r,y) - M{ri,y)M{r2,y)\ 



dy 

(3-9) 
is employed in the more sophisticate model called the MCrcBK model}^ where 
r2 = r — ri . is obtained from the Fourier transform of the numerical results of the 
rcBK evolution equation: 

<PAAk,x,r^)= /;; f d^re-'^-^Vl{l-M{r,y)f (3-10) 

as{k) (271^ J 

where Cp = {N^ - l)/2A^c, y = log(2;o/x) and xq = 0.01, In MC-KLN model, x 
dependence is determined by Eq. p-6p . but in rcBK, x dependence can be obtained 
from the equation. Therefore, we expect that MCrcBK model has more predictive 
power than MC-KLN model. See Eefs. fTTl) . fTTO]) . [19]) . 121)) . fTTT]) for the predictions of 
hadron productions in nuclear collisions using rcBK solutiuons. 

One can use a Gaussian shape for nucleons^^^'"^^^''^^ in the Monte-Carlo 
method. This smooth density profile for a nucleon may be significant for the simula- 
tion of event-by-event viscous hydrodynamics. ^^^' In this case the thickness function 
for a nucleon is given by 

rp(r) = ^exp[-rV(25)]. (3-11) 

The probability of a nucleon-nucleon collision P{b) at an impact parameter b is then 
taken to be 

P(6) = 1 - exp[-kTppib)], Tpp{b) = f dhTp{s)Tp{s - b) . (3-12) 

where (perturbatively) k corresponds to the product of gluon-gluon cross section and 
gluon density squared. We fix k so that integral with respect to the impact parameter 
becomes the nucleon-nucleon inelastic cross section a^N at the given energy: 

cTNNiV^) = I d% (1 - exp[-k{^) Tpp{b)]) . (3-13) 



15 

Note that P{b) broadens with increasing energy, even as the size B of the hard 
valence partons is fixed. 

We finahy note that fluctuations from gluon production can be included in the 
modei^''^ (for MC-Glauber model see Ref. 16911 ) . which is very large contribution 
to the initial higher order anisotropics. 

§4. Energy Loss of Jets inside the Expanding Fluid 

Jet quenching is one of the promising tools to diagnose dense matter created in 
relativistic heavy ion collisions. Energetic partons created in hard scatterings are 
subject to traverse dense medium. Through interaction between these partons and 
medium, various consequences are predicted to take places: suppression of yields for 
high pt hadrons, energy imbalance of a jet pair and so on. In this section we overview 
mainly how to calculate the amount of energy loss in hydrodynamic backgrounds. 
For details about existing energy loss formalisms such as BDMPS-Z,'^^^ GLV,'^^ 
ASW,E3 HjEl and AMV"^ and comparison among them, see e.g. Ref. |I20). 

First attempt to combine energy loss calculation with hydrodynamic simula- 
tions was made in Ref. I12ip . After that, the hydro + jet model based on mini-jet 
production from PYTHIA^^^' and its propagation in full three dimensional ideal hy- 
drodynamic backgrounds^'-^ has been developed and systematic analysis has been 
performed including pT spectra, di-hadron correlation functions, interplay between 
soft and hard components and suppression in forward rapidity regions.'^ 'l^'^^^ 

In most of energy loss formalisms, the amount of energy loss largely depends 
on the inverse of mean free path A~^, and, in turn, on the medium parton density 
p = (o"A)~^ from the kinetic theory, where a is the cross section between an emitted 
gluon and a parton in the medium. For example, jet quenching (transport) parame- 
ter can be g ~ Ai^/A, where fj, is typical transverse momentum transfer suffered from 
the medium. When an ideal gas of massless quarks and gluons is adopted for the 
QGP equation of state, medium parton density is calculated from thermodynamic 
variables, p oc e"^'^ oc T^^^ In this way, one can utilize hydrodynamic outputs 
along a trajectory of a jet to quantify energy loss in relativistic heavy ion collisions. 
It should be noted that the above relation is valid only for the ideal gas equation 
of state. Since the medium is strongly coupled according to hydrodynamic anal- 
yses of collective anisotropic fiow, non-perturbative definition of the energy loss is 
demanded. 1^23 

Initial transverse positions of jets at an impact parameter b are determined 
randomly according to the probability distribution specified by the number of binary 
collision distribution, 

P{r^,b) oc TA{r^ + b/2)TA{r^ - 6/2) (4-1) 

Initial longitudinal position of a parton is approximated by the boost invariant dis- 
tribution: rjs = y, where y = (1/2) ln[{E + pz)/{E — pz)] is the rapidity of a parton. 
Jets are freely propagated up to the initial time tq of hydrodynamic simulations 
by neglecting the possible interactions in the pre-thermalization stages. Jets are 
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assumed to travel with straight line trajectory in a time step: 



^^* = Tk? ^^^' (^ = ^'y)' (4-2) 

rn-j-cosh^y — rjs) 

Ar]s = - tanh(y - r]s)AT, (4-3) 

r 



where rriT = dm? +p^ is a transverse mass. We obtain the total amount of energy 
loss for a sample jet 

dr/>(a;(T),r)(r-ro)lnf^j (4-4) 

when we employ the approximated first order formula in the opacity expansion P^^ 
Here tq is the initial time of hydrodynamic simulations, Eq is the initial energy 
of a jet which is Lorentz boosted by the four flow velocity as PqU^. The formula 
roughly gives L^ dependence of energy loss in the case of a static medium, which is 
manifestation of the Landau- Pomeranchuk-Migdal effect. 

It would be interesting to see a response of the medium to parton energy loss. 
This could be neglected at the RHIC energy where initial energy of jets are not so 
large. On the other hand, jets with a few hundred GeV can be produced at the 
LHC energy. If the lost energy is quickly thermalized, one can solve hydrodynamic 
equations with a source term from energy loss. 

(4-5) 

^o), (4-6) 

(4-7) 

where initial position of a jet is (xq, yo, zq) and the jet is supposed to traverse in 
a straight trajectory in the x direction. These equations can be solved numerically 
using the algorithm mentioned in Sec. | 



d^Tf^^ix) = r{x), 








^"-■-f*(- 


-xo - 


-t)5{y- 


- yo)S{z 


J^ = J'^ = 0, 









§5. Transport Model 

Relativistic transport models have been successfully used to simulate from low to 
high energy nuclear collisions. A great advantage of transport models is that one can 
study space time evolution of system even if it is n ot equilibrate d. Transport models 
include, for example, intra nuclear cascade models ,'^^''^^''^^ Boltzmann-Uehling- 
Uhlenberck (BUU) models,li23 RQMD,lEn},[l3D qgsM,!I2II ARC,E3J ART,'i24i UrQMD,li33 
JAM,i36),i37) HSD,[I3S} and GiBUU.^SS The partonic phase in the early stages of 
heavy ion collisions has been studied by the parton cascade models; VNI,'!^ VNI/BMS,'^^''!^ 
ZPC,iHl MPC,E1 GROMITi-^S) and BAMPS.ESJ Partonic cascade in which only 
the elastic scattering is included does not account for the observed elliptic flow for 
typical gluon-gluon pQCD cross section of agg ~ 3 mb.'^^ Inelastic process such 
as 2 —7- 3 scattering has been shown to be very important for the energy loss of the 
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high energy jets as well as the thermalization of the system. IH^'lHS'IHS'lISo} ,1113 ,1113 
The transport models such as the multi-phase transport model (AMPT)^^ and the 
Paton-Hadron-String Dynamics (PHDS)^^ include the dynamics of both partonic 
as well as hadronic phase. In the following, we shall review mainly the transport 
model JAM.E3 

5.1. Scattering algorithms 

The geometrical interpretation of the cross section for collisions of particles is 
employed by various models, i.e., particles scatter when their closest distance is 
smaller than -y/cr/vr, where o" is a total cross section. The particles are propagated 
along classical trajectories until they scatter or decay. The scattering of two particles 
is determined by the method of closest distance approach; two particles scatter, if 
the impact parameter (closest distance) 6rei for a pair of particles becomes less than 
the interaction range specified by the cross section: 



6,, < f-^ (5-1) 

where o^{^/s) is the total cross section for the pair at the cm. energy ^/s which de- 
pends on the incoming particle pair. This collision method has been widely used to 
simulate high energy nucleus- nucleus collisions. Because of a geometrical interpre- 
tation of the cross section, this method violates causality, and time ordering of the 
collisions in general differs from one frame to another. Those problems have been 
studied by several author^Eli''iS3 

Our approach is motivated by Kodama's workP^'l^^ in which Lorentz invariant 
expressions are obtained by considering the rest frame of one of colliding particles. 
The closest distance 6rei is defined by the distance in their common cm. frame. 
Suppose that the coordinates and momenta of two colliding particles are denoted by 
xi = {ti,xi), X2 = it2,X2), pi = {Ei,pi) and p2 = (£'2,^2)- The trajectories of 
particles are 

xl{t*) = xl{tl) + vl{t*-tl), 

xl{f) = xl{t*2) + vl{f-t*2), (5-2) 

where asterisks represent quantities in the two-body cm. frame. The time of closest 
approach t* may be obtained by the condition that the relative separation becomes 
perpendicular to the relative momentum, (a;^ — X2) ■ {p\ — P2) = 0, 

,. ,*_ v*-{x*-vl{t\-tl)) _ v*-{x\{t\)~x%{t\)) 

,. ,*_ V*-{x* + vl{tl-tl)) V* . jxlit*,) - xUt*2)) ,,. 



where x* = a;|(f|) — 0^2(^2)) t* = ^1 — *2; ^^^ '"* = '"i—'"2- Then the closest distance 
ftrei is expressed as 

^2 ,*2 



{x* ■ v*)"^ 

5rel = X* -Z2 (5-5) 
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One can express 6rei in terms of Lorentz invariant scalars 



ft, = -.= + ^ + ^ (5.6) 

by using the following variables 

X = Xi - X2 = {ti - t2, Xi - X2), (5-7) 

P = pi-P2 = (El- E2,pi-P2), (5-8) 

P = Pi+P2 = iEi + E2,pi+P2), (5-9) 

'i = P-^P^ (5-10) 

where we have used the Lorentz invariant expressions for —{x* •p*)'^ /p*"^ = {q-x)"^ /q^ , 
and transverse square distance — x^ + (P-x)'^/P^ corresponds to the squared distance 
x*"^ in the two-body cm. frame. The following expression for impact parameter in 
terms of the original momenta will be useful 



^2 _ ^2 Pi{P2 ■ xf + pl jpi -x)^ -2{pi ■P 2){pi ■ x){p2 ■ x) 

{pi ■P2Y -P1P2 



"rel — —X — —:^ -3-2 • l^'J-J-J 



One can obtain Lorentz invariant expression for the collision time directly from 
Eq. ()5-4p . but we can start from the covariant equation of motion. By using the 
relations El = {P ■pi)/\P\ and t* = {P ■ Xc)/\P\., Eq. ()5-2p can be rewritten by using 
4- vectors: 



x\ 



Xl + 
X2 + 


P 


■{Xc 


-Xl) 


P 


P- 

■{Xc 


Pi 

Pi 
-X2) 

-P2 



P-P2 

where x'^ and x'2 denote the position of the closest approach. Using the condition 
{x'-^ — x'2) ■ (pi — P2) = 0, one can get the collision times as 



P -{Xc- Xl) _ {P ■ x){p2 ■ p) - {p ■ x){P ■ P2)_ 

y 

(5-12) 



P-Pl {P ■ Pl){P ■ P2) - {P ■ P2){P ■ Pi 

P ■ {Xc- X2) {P ■ X){pi ■ p) - {p ■ X){P ■ Pl) 



P-P2 {P ■ Pl){P ■ P2) - {P ■ P2){P ■ Pl)' 

Two colliding particles are propagated to each collision point of the closest approach: 

pUx - Pi) - (Pi •P2)(a : • P2) , 

(piP2)^ - pIpI 
Pi(x • P2) - (Pi ■P2)ix ■ Pl) 

(piP2)^ - pIpI 



xi{tc) = xi{ti) -] 7Z-r~^ Tl. ?'!' (^'13) 



X2{tc) = X2{t2) 7Z~r\^ 2"^ P2- (5T4) 



The two collision times are in general different, because of the finite spatial separa- 
tion. Therefore one needs a prescription to choose the time of the ordering of each 
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collision. In the model, we assume that the collisions are ordered by the average 
time torder = (^i.coi + *2,coi)/2 in the computational frame. Other choice is to use the 
time of closest approach for the two colliding particles in the computational frame 
as used in UrQMD, 

tcow = 7 ^ (5-15) 

The effect of different definition of ordering time was investigated in Ref . :153p . The 
problem of superluminous signals has been studied in Ref. I155P 

In this geometrical method for two body collision, collision occurs in the sep- 
arated points, which is the same as action-at-distance interaction. In order to 
solve this problem, one can use the "full-ensemble" method or "subdivision" tech- 
nique^^ '^^^'l^^ in which cross section is reduced by the factor of the number of test 
(over sampling) particle a/Ntest in order to recover the local nature of Boltzmann 
collision term in the limit of A'test — ^ oo- However, this method is in general compu- 
tationally expensive for large A'test- A faster method which is called "local-ensemble" 
method has been proposed in Ref. I157p . 

In order to recover the problem of collision ordering, the stochastic niethocpSS|i.ll5l|i,lll6( 
can be used in which probability is used to determine the collision instead of geo- 
metrical interpretation. The collision probability for two particle collision during the 
time interval At in the volume element V is given by 

a At , , 

where v^ei is the relative velocity of the scattering particles. In the limit of Attest — ^ 
cxo, Z\t — )• 0, y — )• 0, the stochastic algorithm will converge to the exact solutions of 
Boltzmann equation, and as a result, recover the Lorentz invariance. Inclusion of 
three body collision in the stochastic method is straightforward. 

5.2. Cross sections 

In this section, we summarize modeling of various hadron-hadron (hh) collisions 
in JAM. the inelastic hh collisions produce resonances at low energies while at high 
energies color strings are formed and they decay into hadrons according to the Lund 
string model with some formation time. Formation point and time are determined 
by assuming yo-yo formation point .1^^ This gives roughly formation time of 1 fm/c. 

5.2.1. Resonance productions 

In our approach, it is necessary to input various hadron-hadron (hh) cross sec- 
tions. We use parameterized cross section for total and elastic collisions. Total 
hadronic cross section in JAM is divided in general by various processes: 

O-tot(s) = O-cl(s) + O-ch(s) + (^a.nn{s) + O-t-R(s) + <Ts-r(s) + <Tt_s(s) + (Ts-s{s), (547) 

where (Tei(s),(Tch(s) and crann(s) denote the elastic, charge exchange and annihilation 
cross sections, respectively. (Tt-R(s) and crs-R(s) are the t and s-channel resonance 
production cross sections, such as NN — )• NA or vrA^ — )• A, and crt-s(-s) and o"s_s(s) 
are the t and s-channel string formation cross sections. 
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Particle production in hh collision is modeled by the resonance formation for 
discrete resonance region, 

hih2 ^ hihl, hih2 ^ h\h*2, or /ii/i2 o h%. (5-18) 

Resonance regions for baryon-baryon (BB), baryon-meson {BM) and meson-meson 
collisions are assumed to be y/s < 4, 3 and 2 GeV, respectively. In JAM, cross 
sections of various resonance formation in BB are parameterized to reproduce pion 
multiplicities. In UrQMD,'^^ BUU,'^® matrix elements are fitted to the available 
data for pion production cross sections. 

The cross section for the inverse process such as h\h2 — )• /11/12 can be obtained by 
the detailed balance formulsP^''^^''^^ which takes the finite width of the resonance 
mass into account. The differential cross section for the reaction (3,4) — )• (1,2) can 
be expressed by the cross section for (1, 2) — )■ (3, 4); 

^^34^12 ^ (2^1 + 1)(252 + I) pf2 ^^12^34 1 

dn (253 + l)(254 + l)p34 dQ J J p^AA{m^)A{mi)d{ml)d{ml) ' 

(5-19) 
where Si denotes the spin of incident or outgoing particles. Mass distribution func- 
tion A{mf) for resonances with the mass rriR is given by the relativistic Breit-Wigner 
function 

A{m ) = —- - 2 \2 I 2~F7 V2' (5-20) 

A/ (m^ — rrij^y + mj^i [my 

where AT denotes the normalization constant. 

The resonance formation cross section for MB and MM collisions is computed 
by using the Breit-Wigner formulsPSSJ-lHOJ (neglecting the interference between res- 
onances). 



a{MB ^R) = ^^^ Yl \C{MB, R)\^ 

R 

i2SR + 1) rR{MB)rRitot) 



"cm 



(25m + 1)(25b + 1) (^/i - mny + r^(tot)2/4 ' 

Sr, Sb and Sj^j denote the spin of the resonance, the decaying baryon and meson 
respectively. The sum runs over resonances, R = A^(1440) - A^(1990), Z\(1232) - 
Z\(1950), yl(1405) - yl(2110), i7(1385) - 17(2030) and ="(1535) - .= (2030). Actual 
values for these parameters are taken from the Particle Data Group^^' and adjusted 
within an experimental error bar to get reasonable fit for MB cross sections. The 
momentum dependent decay width is used for the calculation of the decay width in 
Eq. dEHD, 






Pcmsirnji) J 



where i and Pcms("T') are the relative angular momentum and the relative momentum 
in the exit channel in their rest frame. The Breit-Wigner formula Eq. (|5-21|) is used 
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for meson-meson collision for resonance production as well. Meson resonance states 
are included up to about 1800 MeV. 

When there is no available experimental data for cross section, we use additive 
quark model™ '133 

atot = a.j,^^ (l - OA^) (l - OA^] , (5.23) 

where ctnn is the total nucleon-nucleon cross section, rii is the number of constituent 
quarks in a hadron, and n^j is the number of strange quarks in a hadron. Cross 
sections involving hadrons with many strange quarks such as <j) = (ss) or f2 = (sss) 
are suppressed. Interestingly, this leads to violation of mass ordering in differential 
elliptic flow, e.g., for cj) mesons.l^SSJ This expression is a good approximation above 
the resonance region where cross section becomes flat. Additive quark cross section 
yields Cx-p ~ 21 mb and (Tj\p ~ 35 mb, which are consistent with the experimental 
data. 

5.2.2. String formation 

At an energy range above \/i > 4-5 GeV, the (isolated) resonance picture breaks 
down because width of the resonance becomes wider and the discrete levels get closer. 
The hadronic interactions at the energy range 4-5 < ^/s < 10-100 GeV where it is 
characterized by the small transverse momentum transfer is called "soft process", 
and string phenomenological models are known to describe the data for such soft 
interaction well. The hadron-hadron collision leads to a string like excitation lon- 
gitudinally. In actual description of the soft processes, we employ the prescription 
adopted in the HIJING modeP^ to treat soft excitation processes. In HIJING or 
FRITIOF,^^ excited strings after interaction have the same quark contents unlike 
the model based on Gribov-Regge models such as Dual Parton Models (DPM)ffSz} 
or the VENUS model.Q® In DPM or VENUS, strings typically have different quark 
content than original hadrons as a result of color exchange. 

We shall review string excitation employed by HIJING. In the center of mass 
frame of two colliding hadrons, we introduce light-cone momenta p = E ± pz- 
Assuming that beam hadron 1 moves in the positive z-direction and target hadron 
2 moves negative z-direction, the initial momenta of the both hadrons are 

Pi = (p^pF^Ot), P2 = {pt,P2^0T) ■ (5-24) 

After exchanging the light-cone momentum {q^,q^,PT), the momenta will change 
to 

Pi = (Pi +P2 -P}^ P/' Pt), P2 = (P/> P\ +P2 -P'}^ -Pt), (5-25) 

where final momenta pt and p7 are related to the longitudinal momentum transfer 
q^ as 

p^=q^+pt, P]=q'+Pi- (5-26) 

Using light cone momentum transfer x^ defined by 

x+ = ^, x- = ^, (5-27) 
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Final momenta are given by 

Pi = ((1 - x^)Vs, X" Vs, pt), P2 = {x'^Vs, (1 - x~)y/s, -pr). (5-28) 
Thus the string masses wih be 

Mf = x-{1-x+)s-pI, Ml =x+{1-x~)s-pI, (5-29) 

respectively. Minimum momentum fractions are x^j^ = j?^/P^ and x~j^ = p^ /P~ . 
For the probability for light-cone momentum transfer in the non-diffr active events, 
we use the same distribution as that in HIJINGP^ 

^(^ > - (^±2 + ,2/,)i/4 (5 30) 

for baryons and 

^^"^ ^ " (2;±2 + c2/s)i/4((i - x±)2 + cysy/^ ^^'^^^ 

for mesons, where c = O.lGeV is a cutoff. For single-diffractive events, in order to 
reproduce experimentally observed mass distribution dlVP /M'^, we use the distribu- 
tion 

The same functional form as the HIJING modeP^ for the soft px transfer at 
low pt < Po is used 

/(Pt) = {(p?^ + c?)(p|.+pg)(l + e(P^-P°)/^^)}"' , (5-33) 

where ci = 0.1 GeV/c, po = 1.4 GeV/c and C2 = 0.4 GeV/c, to reproduce the high 
momentum tail of the particles at energies -Eiab = 10-20 GeV. 

5.2.3. String decay 

The strings are assumed to hadronize via quark-antiquark creation using Lund 
fragmentation model P YTHI A6 . 1 \^^ Hadron formation points from a string frag- 
mentation are assumed to be given by the yo-yo formation poinlP^ which is defined 
by the first meeting point of created quarks. Yo-yo formation time is about 1 fm/c 
assuming the string tension k = 1 GeV/fm. 

In the Lund string model, space-time coordinates and energy-momentum co- 
ordinates for the quarks are directly related via the string tension.1^''i^ Let us 
consider one-dimensional massless gg string in the cm. frame. If x^ = ti^Xzi denote 
the light-cone coordinates of the i-th production point, then the light cone momenta 
Pj^ = En ± pzi of the i-th rank hadron which is produced by the energy-momentum 
fraction Zi from (z — l)-th string p^ = Zipf_^ are fixed by 
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constituent formation point 
yo-yo formation pointy 




Fig. 2. Space-time picture of the qq string motion and the definition of formation points. 

with initial value for quark qq moving to the right 

+ W 

where W corresponds to the string initial invariant mass. The probability distribu- 
tion for the momentum fraction Zj is given by the Lund symmetric fragmentation 
function 



f{z)^^—^e^v[-h- j^^l , (5-36) 

where a and b are parameters which have to be fitted to experimental data, and 
m and p± denote the mass and transverse momentum of the produced hadron, 
respectively. This form can be obtained by the condition of the left-right symmetry 
for the decay of string. 

Using the relation pfp~ = 'mf_^ with 7raj_|_ being transverse mass of the i-th 
hadron, we have the recursion formulaelSSJ 

mjj \ 2 1 - Zi 1 



a^i =(l-^i)a;«-i, X. =x._i + (-— ) — ^. (5-37) 

Yo-yo formation points where the two qq meet for the first time are obtained as 

xr'° = (xti,x-), (5-38) 
and constituent formation points where qq is created 

:^^ = {xt,xr). (5-39) 

In RQMD,'i21J' the formation points of hadrons are calculated as the average of the 
two qq production point as 

^RQMD ^ ht+^U ^ ^^^1 . (5-40) 
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Clearly, one can see that 

It is assumed in UrQMD and JAM that the yo-yo formation point is assigned to 
the space-time point for produced hadrons. The produced hadrons which have the 
original constituent quarks are allowed to scatter with reduced cross section accord- 
ing to the number of constituent quarks. For example, if baryon has one original 
constituent quark, cross section during the formation time would be reduced by the 
factor one-third. 

§6. Non-Abelian plasmas 

Understanding the dynamics of the non-equilibrium system and the process to- 
ward the thermalization is a very important outstanding question in high energy 
heavy ion collisions. In Ref. 1171]) . "Bottom- up thermalization" was proposed based 
on the Boltzmann equation with inelastic processess. Later, full numerical simula- 
tions of the Boltzmann equatiorP^ show that inelastic gluon scatterings are very 
important processes for thermalization of gluonic systems. It was, however, pointed 
out that the soft classical color fields may play an important role in the dynamics 
of thermalization. Specifically, Weibel-like QCD plasma instabilities may develop 
due to anisotropic distribution of hard plasma particles (partons).'^''^''^ First 
numerical simulation was done within the hard-loop approximation by solving lin- 
earized Vlasov equation in one-dimension.^^ In this analysis, exponential growth of 
gauge fields due to Abelianization of the field was found. In full 3-dimensional sim- 
ulations,'^^ 'l^^ it was found that instabilities grow linearly once the field strength 
becomes large and non-Abelian self-interactions among gluon fields become non- 
perturbative in the case of moderate anisotropic momentum distribution. In this 
case, non-Abelian interactions develop a cascade of energy from soft modes to hards 
modes.'^^ However, in the case of extreme anisotropy,'^^''^^ energy coming from 
plasma particles due to Weibel-like instability does not go to the soft mode, instead 
go back to the hard scale very rapidly which is called "ultraviolet avalanche" in 
non-Abelian plasmas. 

Instabilities of the gluo n fields have bee n also investigate d by employing pure 
classical Yang-Mihs equation.Il^'E3>[IHl.lIHD it was pointed outPS' that instabilities 
in the classical Yang-Mills field are the Nielsen-Olsen type instability. Numerical 
analysis was performed in Ref. I182p . In Nielsen-Olsen instability, unstable mode 
exists at zero momentum, while zero momentum mode is stable in Weibel instability. 
This is one of the crucial differences between Weibel and Nielsen-Olsen instabilities. 

In this section, we review how to solve Wong- Yang-Mills equation based on the 
technique developed in Ref. I177|) below. 

6.1. Simulation of non-Abelian plasma 

We present numerical method of the classical Vlasov transport equation for 
gluons with non-Abelian color charge q*^ 

pf'id^ + gq'^F^.d^p + gfatcAlq'd,a]f{x,p, q) = (6-1) 
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where f{x,p,q) denotes the single-particle distribution function and g is the gauge 
coupling. The Vlasov equation is coupled self-consistently to the Yang-Mills equation 

D^F^^'^ = r =gl^dqqv''f{t,x,p,q), (6-2) 

with v" = (l,p/p), D^ = d^ + igA^ F^^ = df,A^-duA^ + ig[A^,A^], and J" denotes 
the current generated by plasma p arti cles. This set of equations reproduces the 
"hard thermal loop" effective theor}U^ near equilibrium. 

The Vlasov equation can be solved by replacing the one-particle distribution 
function with many "test particles" : 

fix,p, q) = ^ Yl ^(^ - ^^(tMP - Piii^Q - 'l^ii))^ (6-3) 

jVtest 

where iVtest is the number of test particles, and Xi(t), Pi{t) and qi{t) are the coor- 
dinate, momentum and the color charge of the test particle. Prom this assumption, 
one gets the Wong equatiorP^ for the i-th "test particle" 

^=v., ^=gqt{E'^ + v,xB^), ^ = -igvn\,q^] (6-4) 

The time evolution of the Yang-Mills field can be followed by the Hamiltonian 
methocP^ in the temporal gauge. The Hamilton equation of motion for gauge field 
Ai = A^t"" and the color electric field Ei = Eff^ under the temporal gauge ^'^ = is 

d./^'i J dill V ^ „ ^A , •■. , : 

-^=E^ -^ = Y.^3^i'-'^^ (i,i = x,y,z) (6-5) 

i 

One may descritize the equation on the lattice size a by introducing the link 
variable Ui{x) = eyip{iagAi{x)) 

U,{x)=iagE\x)Ui{x), (6-6) 

Ki^) = ^E^r {i''[f^..W - U]{x-3)U,,{x-j)U,{x-j)]] - J: (6-7) 

where, the plaquette is defined as Uu = Uij{x) = Ui{x)Uj{x + i)Ul{x + j)U:{x), and 
t"- = a"" is the Pauli matrix for SU(2) case. This equation is covariant under the 
lattice gauge transformation 

Ui{x) -^ V{x)Ui{x)V\x + i), E\x) -^ V{x)E'{x)V^{x), (6-8) 

Typically one uses the time step size of At ~ 0.05a to ensure energy conservation 

i,x n j 

and Gauss law 

Y, (E^{x) -UJix- i)Ei{x - i)Ui{x - i)) = p{x) (6-10) 
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Numerical method has been developed to solve the equations ()6-4p in Ref. I186p 
which is the non-Abelian extension of the nearest-grid-point (NGP) method. In the 
NGP method, the charge density is obtained by counting the number of particles 
within a cell. A current J{x) = Q5(t — Across )/-^tcst is generated only when a particle 
with the color charge Q crosses a cell boundary from x to x -|- i at the time tcross- In 
order to satisfy the requirement of the lattice covariant continuity equation 

p(x) = Y^ul{x- i)Ji{x - i)Ui{x -i)- Ji{x) , (6-11) 

i 

the color charge must be parallel transported 

Q(x + i) = ul{x)Q{x)Ui{x) . (6-12) 

At each time step, an effect of magnetic field, which causes a rotation of momentum 
and does not change the energy of the particle, on the particle is taken into account. 
But the magnitude of the momentum can change when it crosses a cell boundary. The 
final momentum |pfin| after crossing the cell can be fixed by the energy conservation 

I I , -'*tcst ,-,2 I I , -' * test / „ T\2 /n i oN 

|Pmi| + ^— £;;„; = IpfinI + ^—(-Bini- J) • (6-13) 

Substituting J = Q/A'test and when a particle crosses the cell boundary in x- 
direction, one obtains for new momentum in the x-direction 



Px = V IPinil - ^xsiniQ + QV(2A^test) " pI , (6-14) 



where Ipdnl = \/Px~^P'± is used. As we explained above, because of the current 
discontinuities introduced in NGP method, large amounts of noise are generated. In 
order to eliminate this noise, we need a large number of test particle. NGP method 
was applied ijJlaDiU^'llSSi for one-dimensional case in which gauge field is assumed 
to depend only on one-direction: A(x,y,z) = A(z). We call this approximation 
as ld-3v simulation because particle velocity has three direction. In ld-3v case, we 
can take a very large number of test particle to obtain sufficiently smooth current. 
However, for full 3d-3v simulation, computational cost may be very expensive due 
to the larger number of test particle required, and we need improved algorithms. 

6.1.1. Particle-in-Cell (PIC) simulations 

In Abelian plasmas, particle-in-cell (PIC) method is widely used in which smoothed 
currents are employed^^*^''^^^''^^^' in order to suppress numerical noise. The rn--th 
order shape function (b-spline) Sm can be obtained by the convolution of 5m- 1 with 
the nearest-grid-point weighting function Sq 

/oo 
So{x')Sm-i{x - x)dx' (6-15) 

where zero-th-order function is the flat-top function 

5o(x) = | ^ fo^l^l^2' (6-16) 

I otherwise 
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Commonly used is the first-order function Si wliich corresponds to linear interpola- 
tion 

1 — |x| for |a;| < 1, 



\ otherwise 

The charge density at each lattice site is given by the superposition of a particle, 

^("■•') = s, 4iizi^- (^) '- (^) «- (^) ■ <"«' 

where (x, y, z) is the coordinate of particle, and Ax, Ay and Az are the lattice spacing 
in the x, y and z directions, respectively. However, discrete continuity equation would 
not necessarily be satisfied on the lattice. One may need to solve Poisson equation 
to recover Gauss law. More efficient numerical methods that satisfy exactly the 
lattice continuity equation have been develop ed \^^ ^^ 'l^^ In PIC simulation, the 
total amount of charge is distributed over its surface, and each charge contributes 
to the charge density in several cells (8 cells for first order in 3-dimension) . In NGP 
method, current is generated only when a particle crosses the cell boundary, but in 
PIC, current will be continuously generated in a given time interval which is the 
amount of charge crossing a cell boundary. 

Let us try to construct a current which satisfies the lattice continuity equa- 
tion. Difference of charge density within a time step At at a lattice site (i,j) in 
2-dimensional case is 

P*+^*(^,j)-p*(^,i) 

= E ^^{ST{<t + ^t))ST{y{t + At))-ST{x{t))ST{y{t))} 

particle 

ft+At 1 
= E y j^{qST{x(t))ST{y{t)))dt, (6-19) 

particle 

where SY^{x) = -^Sm (^^)- By using the identity for the derivative of Sm 

"'^"Ir/^'^^ = -^[Sn.-iix + Ax/2) - 5„_i(x - Ax/2)] (6-20) 

which follows from the Eq. (|6-15p . one finds that the current defined as 

T /. • ,N 1 v^ r, fx — i — Ax/2\ ^^^ fy — j\ TT, (z — k 

^"(^'^'^) = A^Z^ E '^-"^-^ ( Ax r- \M) ^- [-AT 

pcirticlG 

(6-21) 
satisfies the lattice continuity equation, where 

/t+At 
S„,{x)dt . (6-22) 

We consider a particle moving from {xi,yi) to (x2,y2) during a time step At, 
i.e., X2 = xi + VxAt, 2/2 = 2/1+ VyAt. We also restrict ourselves to the case in 
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which (x2,y2) belongs to the same lattice site {i,j)- In this case, The explicit form 
for the first shape-factor corresponding linear smearing case is given by 

Pihj) = . . (1 - x){l - y), p{i,j + 1) = x{l - y) , 

AxAy AxAy 






J^i^^j) = ^Tzn^M'^ - Wy)^ Mi, J + 1) = -^-^F^Wy , 



•^^^^'^■^ = A^^'^^ " ^^■)' ^^^' + ^'^'^ = A^^y^^ ' ^^-^^^ 



Fx = q :^r- ' Fy = q ^^ . (6-25) 



Wx = ^ « , VFy = — ^ J . (6-26) 



where x = (xi — i)/Z\x, y = {yi — j)/Ay, and {Fx,Fy) = F represents the charge 
flux 

At ' y~'^ At 

Wx/y is defined at the midpoint between the starting point (xi, yi) and the end point 

(2^2,^2) 

2;i + X2 . „. yi + y2 

We use the 'Zigzag scheme' developed in Ref. I195P in case particle crosses the cell 
boundary within a time step. 

Finally, we note that the electromagnetic forces should be smeared in a similar 
way when a particle momentum is up dated .l^^ Consider the time derivative of the 
total energy: 

-^ = -Y^E.J+ Y. q,E{x,{t))-Vi. (6-27) 

lattice particles 

It is clear that the interpolation function for E should be the same as that for J in 
order to achieve good energy conservation in the simulation. The electric field E{x) 
at the particle position x is then obtained from 

E^ix) = Y, SlSlS]E^{i,j,h) , (6-28) 

lattice 

while the magnetic field is given by 

B^{x) = Y^ SiSp?^B^{i,j,k) . (6-29) 

lattice 

This is motivated by the relations E = —VA^ — A and B = V x A. 

Next we consider the discretization of the equation Eq. (|6-4|) for momentum 
update. The difference form of the equation is 

p{t + At/2) - Pit - At/2) , / Pit + At/2) + Pit -At/2) B^it) 

— gq [ E [t) -\ X 



At ^-i y y J ' 2 eit) 

(6-30) 
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Equation ()6-30p can be solved by the Buneman-Boris niethocF^''^^ as follows: 

Pit) = Pit- At/2) + ^Eit), (6-31) 

eit) = \pit)l (6-32) 

. ^ ^t , , Bit) 

p'it)=pit) + —pit)x-^, (6-33) 

^^(^) = ^(^) + l + iB/eit)At/2)^ TP'^'^ ^ ^' ^'-''^ 

Pit + At/2) =p2it) + ^Eit), (6-35) 

where E = gq°'E"' and B = gq'^B"'. This scheme is time reversible and the overall 
momentum integration is accurate to second-order in the time step. 

6.2. PIC simulations in non-Abelian gauge theories (CPIC) 

An extension of the charge conserved method to the smearing method in the 
non-Abelian case has been proposed in Ref. I177P in which the current is defined as 



^^(l-iy^), Jxil,] + l) = Qy^^ 



Ui,j) = Q^^il-Wy), J^iij + l) = Qy^^Wy, (6-36) 



Jyii,j) = Q^^^^il-W,), Jyii + l^j)=Qj-l^W,, (6-37) 

where we define the parallel transport of the charge in two dimension as 

Q. = Ulii,3)QU.,ii,3), Qy = Ulii,j)QUyii,j) . (6-38) 

One can easily check that this satisfies the lattice covariant continuity equation, 

KhJ) = XI ^^(^ ~ x)J^ii - x)Uxii -x) - JxiiJ), (6-39) 

X 

for sites (i, j), (z + 1, j), (i, j + 1): 

pii,j) = .Ui,j) + Jyii,j), (6-40) 

pii + 1, j) = Ulii,j)J,ii,j)U,ii,j) - Jyii + 1, j), (6-41) 

p(i, j + 1) = UliiJ)Jyii,j)Uyii,j) - MiJ + 1), (6-42) 

Equations (|6-40p . ()6-4ip and (|6-42p are consistent with the following definitions of 
the charge densities: 

pii,j)=Qil-x)il-y), (6-43) 

pii,j + l) = Qyil-x)y , (6-44) 

pii + IJ) = Q^xil - y) . (6-45) 

However, since a particle's color charge depends on its path, so does pii + 1, j + 1) 
and we are not able to calculate it from the charge distribution itself. Rather, we 



30 Tetsufumi Hirano, Yasushi Nara 

directly employ covariant current conservation to determine the increment of charge 
at site {i + 1, j + 1) within the time-step. In this way, we can satisfy Gauss's law in 
the non-Abelian case. 

Finally, we have to check that Tr((5^) is conserved by this smearing method. 
This is true when the lattice spacing a is small, as the total charge of a particle is 
given by 

Qo = Q{l-x){l-y) + Qxx{l-y) + Qy{l-x)y+[apQxy + {l-ap)Qya]xy , (6-46) 

where the Op depend on the path of a particle and Qxy = Ul{i,j + l)QyUx{i,j + l), 
Qyx = Uy{i + l,j)QxUy{i + 1, j). If we require that Ti{Qq) be constant, then the 
cross terms, for example Tr{QQx), have to vanish. This is true when a is small, 
because Tr((5[^,Q]) =0: 

Tr{QQx) = Tr(Q(Q + iga[Ax,Q] + ©(a'))) = TV(Q2) + ^(a'). (6-47) 

Therefore, the magnitude of the color charges is conserved for small lattice size a in 
our method. 

§7. Summary 

We have reviewed recent progress on the development of dynamical models in 
heavy ion collisions. Our emphasis was put on the technical aspects of the models. 
Initial particle production processes are obtained by the CGC framework based on 
the kT factrization formula. Vlasov simulation of non-Abelian plasma can be followed 
by the non-Abelian extension of particle-in-cell method. In the locally thermailzed 
stage, relativistic hydrodynamics can be used to describe space-time evolution of 
matter. Energy loss of energetic partons in the expanding fluid elements is necessary 
for the consistent description of the jet quenching as well as two-particle correla- 
tions. Finally, effects of final state interactions during the hadronic gas state are 
also important for the realistic simulation of heavy ion collisions. 

Recently there have been made significant progresses for theoretical approaches. 
Viscous hydrodynamic simulations have been performed by many groups^SS' to ex- 
tract transport coefficients such as a ratio of shear viscosity to entropy density. 

Numerical solutions of the classical Yang-Mills equations are first employed for 
the initial condition of event-by-event hydrodynamical simulations.'^ JIMWLK 
renormalization group evolution equation was used as initial conditions for the clas- 
sical Yang-Mills equation^^^' in order to incorporate the rapidity evolution of the 
probability distribution for Wilson lines. This approach may give an important initial 
condition for hydrodynamic simulations. SU(2) plasma simulations of Boltzmann- 
Vlasov equation including both collision term for hard particles and soft interaction 
by classical Yang-Mills fields are performed in Refs. fT98ll . fT99]) . [2nn|) . [2nT]l . It was 
demonstrated that results are independent of the choice of separation scale between 
hard and soft modes. Inclusion of inelastic process is important for the jet energy 
loss in the early stages of heavy ion collision before thermalization. 

One of the outstanding problems in high energy collisions at RHIC and LHC is 
the missing understanding of non-equilibrium dynamics in the early stages. Recent 
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progress on the investigations how non-Abehan plasmas or Glasma approach equi- 
librium state in heavy ion collisions can be found, e.g., in Refs. [202|) .[203) . 204) . 1205]) . 
Non-equilibrium simulations for these approaches will provide insight into the mech- 
anism of thermalization in heavy ion collisions. 
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